%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Numerische Blochgleichungen mit Relaxation, n = (n11 n22 n12 n21).
%  Photonenecho, Doppler-freie Spektroskopie.
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function Bloch
    figure(1);						clf;
    col='rgbmc';

    t = .1;
    g =.03;
    O =  2;
    
    % equação de Schrödinger
    
    for k=1:5
        psi=[1 0]';
        d =-.1*(k-3);
        o =O;
        w =pi/2/O;

        tt=0;
        for tt=0:t:w+5+2*w+5+w
            o = 2;
            if tt<w+5+2*w+5
                o = 0;
            end
            if tt<w+5+2*w
                o = O;
            end
            if tt<w+5
                o = 0;
            end
            if tt<w
                o = O;
            end

            ns(1)=-sqrt(2)*imag(psi(1)*psi(2)');
            ns(2)=-sqrt(2)*real(psi(1)*psi(2)');
            ns(3)=(psi(2)*psi(2)'-psi(1)*psi(1)');
            psi = Schrd2(psi,t,o,d);
            nt(1)=-sqrt(2)*imag(psi(1)*psi(2)');
            nt(2)=-sqrt(2)*real(psi(1)*psi(2)');
            nt(3)=(psi(2)*psi(2)'-psi(1)*psi(1)');

            subplot(221),plot3([ns(1) nt(1)],[ns(2) nt(2)],[ns(3) nt(3)],col(k),'MarkerSize',10,'LineWidth',1.2)
            set(gca,'xLim',[-1 1],'yLim',[-1 1],'zLim',[-1 1],'FontSize',16,'LineWidth',1.2)
            xlabel('\fontsize{16}Im {\it\rho}_{12}')
            ylabel('\fontsize{16}Re {\it\rho}_{12}')
            zlabel('\fontsize{16}2{\it\rho}_{22}-1')
            text(1,-1.7,1.3,'\fontsize{16}(a)')
            view([130 20])
            hold on
            pause(.01);
        end
    end
    hold off

% equações de Bloch

    for k=1:5
       n =[1 0 0 0]';              tt=0;
       d =-.1*(k-3);		  		o =2;
       w =pi/2/o;
       for tt=0:t:w+5+2*w+5+w
            o = O;
            if tt<w+5+2*w+5
                o = 0;
            end
            if tt<w+5+2*w
                o = O;
            end
            if tt<w+5
                o = 0;
            end
            if tt<w
                o = O;
            end

            ns(1)=sqrt(2)*imag(n(3));
            ns(2)=sqrt(2)*real(n(3));
            ns(3)=real(n(2)-n(1));
            n = Bloch2(n,t,g,o,d);
            nt(1)=sqrt(2)*imag(n(3));
            nt(2)=sqrt(2)*real(n(3));
            nt(3)=real(n(2)-n(1));

            subplot(222),plot3([ns(1) nt(1)],[ns(2) nt(2)],[ns(3) nt(3)],col(k),'MarkerSize',10,'LineWidth',1.2)
            set(gca,'xLim',[-1 1],'yLim',[-1 1],'zLim',[-1 1],'FontSize',16,'LineWidth',1.2)
            xlabel('\fontsize{16}Im {\it\rho}_{12}')
            ylabel('\fontsize{16}Re {\it\rho}_{12}')
            zlabel('\fontsize{16}2{\it\rho}_{22}-1')
            text(1,-1.7,1.3,'\fontsize{16}(b)')
            view([130 20])
            hold on
            pause(.01);
        end
    end
    hold off
end

function psi = Schrd2(psi0,t,o,d)
    H =[0    o/2 
        o/2   d ];
    psi=expm(1i*H*t)*psi0;
end

function n = Bloch2(n0,t,g,o,d)
    I =     1i/2;

    M=[ 0      g     I*o   -I*o 
        0     -g	-I*o    I*o 
        I*o  -I*o -g/2-1i*d  0   
       -I*o   I*o     0   -g/2+1i*d];
    n=expm(M*t)*n0;
end

